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I. INTRODUCTION 

The dynamics of crack propagation is an important 
and long standing mystery in solid-state physics and ma- 
terials science jl] [2j , and in the recent years the physics 
community experienced a rebirth of interest in the prob- 
lem of dynamic fracture. The fundamental basis of to- 
day's understanding of the phenomenon fracture traces 
back to Griffith [3] , who realized that the growth of cracks 
is determined by a competition of a release of elastic 
energy and a simultaneous increase of the surface en- 
ergy due to the advancing crack. The uniform motion 
of a crack is relatively well understood in the framework 
of continuum theories f3HS]- Here, the conventional ap- 
proach is to treat the crack as a front or interface separat- 
ing broken and unbroken regions of the material; propa- 
gation is governed by the balance of the elastic forces in 
the materials and cohesive stresses near the crack tip [3- 

. Many characteristic features of crack propagation are 
nowadays well established by experimental studies [TOl - 
[T7] . As soon as the flux of energy to the crack tip exceeds 
a critical value, the crack becomes unstable and starts to 
branch while emitting sound waves. These phenomena 
are consistent with the continuum theory of sharp crack 
tips, but it fails to describe them, because the details of 
crack growth, in particular the chosen crack path and 
velocity, depend on details of cohesion on microscopic 
scales [18]. Nevertheless, empirical energy balances and 
simple propagation laws that are frequently used in en- 
gineering applications, cannot account for the richness 
of actual fracture phenomena. In particular, they can- 
not predict dynamical instabilities of fast moving cracks. 
The fundamental mechanisms of these instabilities have 
been extremely difficult to elucidate because they appear 
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to result from a non-trivial coupling between dynamical 
phenomena inside the crack tip region, known as process 
zone, and (linear) elasticity, with no clear separation of 
scale between atoms and the system size. 

Large scale molecular dynamics (MD) simulations with 
about 10^ atoms allowed to get deeper insights into the 
growth behavior of cracks Although limited to 

submicron samples and very short timescales, these sim- 
ulations were able to reproduce key features of crack 
propagation like the initial acceleration and the onset of 
instabilities. Nevertheless, a detailed understanding of 
the complex physics of crack propagation, in particular 
aspects of the pattern formation process, still remain a 
major challenge [23 . 

At this level, continuum descriptions, in particular 
phase field methods that avoid dynamical artifacts which 
are associated with the breaking of translational and ro- 
tational symmetry [211 [25], offer a useful and comple- 
mentary perspective on crack propagation as a pattern 
formation process. The past years have seen intense ac- 
tivities in phase field modeling of crack propagation (see 
[26] for a recent review) and of defects in general [27] . 

Here, we propose a continuum description of crack 
propagation in the spirit of interfacial pattern formation 
processes. Inspired by the discovery of the Asaro-Tiller- 
Grinfeld (ATG) instability [^HHSO] . we understand frac- 
ture as late and highly nonlinear stage of this elastically 
driven interfacial instability. In its early stage, a linear 
stability analysis of a solid surface under uniaxial load 
reveals that long wave morphological perturbations are 
unstable in the sense that they lead to a decrease of the 
total free energy. Finally, in a later stage of the instabil- 
ity one observes the formation of deep notches, which are 
similar to cracks (see e.g. [3TH53])- Nevertheless, if solely 
accounting for linear elasticity, this instability leads to 
a breakdown of the physical description due to a finite 
time cusp singularity: After a finite time, the unstable 
deep grooves advance with infinitely high velocities and 
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vanishing tip radii (see e.g. [M])- 

This problem can be solved, for instance, by the inclu- 
sion of elastodynamic effects which restore the selection 
of the steady state tip radius and velocity. Based on this 
recognition, a minimal continuum theory of fracture was 
developed using only well-established thermodynamical 
concepts j35j. In this picture, a full modeling of a prop- 
agating crack not only determines the crack speed but 
also the entire crack shape and scale self-consistently, 
which leads to a description as a moving boundary prob- 
lem. The latter was then solved by basically two differ- 
ent methods: First, a sharp interface multipole expansion 
technique |36j and the fully dynamic phase field method 
[37l l38] . Remarkably, already this single parameter min- 
imum model selects steady state propagation velocities 
appreciably below the Rayleigh speed and shows a tip 
splitting instability for high applied tensions. A short- 
coming of the model is a decaying velocity as function 
of the driving force over a significant range of applied 
stresses. 

Recently, a similar continuum model of fracture was 
proposed in |39j , curing the problem of the finite time sin- 
gularity by viscoelastic damping. Apart from the usual 
dissipation directly at the crack surface, viscous bulk 
dissipation takes place in an extended zone around the 
crack. Hence, the incoming flow of elastic energy is par- 
tially converted to surface energy, in order to advance the 
crack, and thermal energy due to viscous damping. How- 
ever, this model does not capture a branching instability 
for high applied loads in case of mode I loading. 

The purpose of the present paper is twofold: First, it 
summarizes and extends the aforementioned work for the 
limiting cases of pure elastodynamics and viscoelasticity. 
Second, it introduces a description which contains both 
effects, therefore capturing the benefits of them and over- 
coming the limitations. We apply the model mainly to 
mode I, but also mixed mode loadings consisting both of 
mode I and mode HI. Two different material transport 
mechanisms are considered and compared. 

The paper is organized as follows: In section [H] we 
present the continuum model of crack propagation in 
elastodynamic and viscoelastic media. Then, in section 
|III[ the crack tip selection principles are discussed. The 
arising free boundary problems are solved numerically by 
the use of sharp interface and phase field methods, as pre- 
sented in section [TV] Finally, we discuss the predictions 
of the model in section Ivl 



II. CONTINUUM MODEL OF FRACTURE 
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FIG. 1: Sketch of crack propagation mechanisms, (a) Bond 
breaking picture, depicted in the reference frame. The black 
line indicates the crack, and there mechanical boundary con- 
ditions have to be applied. The atoms which change their 
atomic configuration due to bond breaking are shown in red. 
(b) The same crack in the deformed configuration, showing 
the separation of interface atoms, (c) Dislocation emission 
from the crack tip leads to blunting. The atomic neighbor- 
hood relations change in the bulk, as illustrated by the green 
atoms, (d) Sketch of the atomic configuration in the reference 
frame for a crack with finite tip radius, (e) The same in the 
deformed state, (f) Propagation of a crack with finite tip ra- 
dius demands mass transport, since atoms have to be removed 
from the tip region. Here we illustrate the crack growth by 
surface difi^usion. The transparent background shows the con- 
figuration in the previous timestep, the solid red atoms the 
interface atoms after advance by one lattice unit. 



We propose a description of fracture in the spirit of 
elastically driven interfacial pattern formation processes. 
In contrast to classical descriptions, where the tip is 
treated as a singular point followed by a mathematical 
cut, we assume the crack to be macroscopically extended, 
and, even more important, to have a finite tip radius 
ro ^ h, a.s can be seen in Figs, [l] [2] and|3] 



This finite tip size implies, that the volume inside the 
crack is also finite, and - depending on the growth mecha- 
nism of the crack - also a description of this inner "phase" 
is necessary. The shape of the crack itself is not an input 
to the model, but is determined self-consistently by the 
equations of motion for the entire crack. In this sense. 
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the description differs substantially from classical mod- 
els, where only equations of motion for the singular crack 
tip have to be postulated or derived. One of the advan- 
tages of such a description is that the entire crack shape 
is a degree of freedom for the model, and therefore not 
only the advance of the crack itself is described, but also 
deformations of the crack contour behind the tip, and - 
what is much more important - path selection is auto- 
matically contained in such models 40 . 

The equations of motion for the crack depend on the 
local elastic deformations, but also the local curvature 
of the interface. In this way, they naturally capture the 
effect of stress release through crack propagation, but 
also the increase of interfacial energy due to crack elon- 
gation, which is the basis for the Grifhth criterion. It 
turns out that the desired self-consistent selection of the 
crack shape is a nontrivial step, since by the aforemen- 
tioned ATG instability the tip tends to become sharper 
and sharper, without a self-consistent selection of a tip 
scale, if only static linear elasticity and interfacial energy 
are taken into account. 

Another important aspect is related to the definition 
of the crack shape. We describe all patterns here in the 
Lagrangian reference frame, which means that in this 
configuration the mechanical deformations are not taken 
into account. In this reference frame a straight cut, which 
is frequently used for a mathematical description of a 
crack, would just appear as a line, and has a vanishing 
tip radius. Under deformation, however, the crack sur- 
faces separate, in particular the distance between the lips 
would scale as Am ~ Kir^^'^/E, where r is the distance 
from the crack tip, Kj the mode I stress intensity factor 
and E is the elastic modulus of the material. Pure elas- 
ticity describes only the deformation of materials, but 
it does not provide evolution equations of the motion of 
the crack. In particular, linear elasticity would predict 
a tr ^ Kjr^^/'^ stress singularity at an infinitely sharp 
crack tip. Physically, one would expect regularization of 
this singularity either by nonlinear effects or a finite tip 
radius ro, which serves as a cutoff parameter. In this 
work we do not follow the first regularization approach 
and consider only linear mechanical models, which will 
be linear elastodynamics and a linear (Kelvin) viscoelas- 
tic model; instead, we consider situations with a finite 
tip radius. 

Let us briefly contrast this description with conven- 
tional pictures of crack growth on an atomistic level. 
In brittle materials, the intuitive interpretation of crack 
propagation is due to the breaking of bonds at sharp 
crack tips. This changes the neighborhood configuration 
for an atom at the crack tip only in the sense that con- 
nections to some of the adjacent atoms is lost, and simul- 
taneously the distances to the other atoms are changed 
due to elastic deformations. This breaking of bonds cor- 
responds to the advance of the sharp cut in a continuum 
model (see Fig. [T^ and b). Ductile effects leads to the 
emission of dislocations from the crack tip, and lead to 
blunting (see Fig. fTl;). These plastic events introduce 



also configurational changes in the bulk in the sense that 
the neighborhood relations are changed. We point out, 
that this is a bulk effect, which of course also reaches 
the crack, since the dislocation lines have to terminate at 
surfaces. Notice that on the continuum level such model- 
ing requires either equations of motions for dislocations 
or - in a coarse grained framework - plasticity models. 
In this article, we will focus on yet another effect, which 
is not captured by the above picture, which is due to 
material rearrangement at surfaces (see Fig. [ijl-f). It 
means, literally speaking, that atoms are removed and 
attached to the crack surfaces at different places, and 
therefore the neighborhood configurations are changed 
now in the sense of a surface effect. For a crack that 
has a finite tip radius already in the reference frame, ad- 
vance of the crack means that atoms have to be removed 
from the tip. They can be deposited again on different 
regions of the crack surface, and in this case we can in- 
terpret the material transport as a surface diffusion pro- 
cess. Alternatively, we could imagine that the removed 
atoms become part of a "gas phase" inside the crack. A 
gas has of course a lower density than the solid, which 
would require that ultimately the gas atoms have to be 
ejected from the crack (if the density difference is not 
accommodated by the crack opening). For convenience, 
we will not consider this fast "hydrodynamic" transport 
and ignore the density difference between the solid and 
the "dense gas phase" . Notice that in both cases of ma- 
terial transport the atoms undergo long-range transport 
(on the atomic scale), and therefore their neighborhood 
configuration changes completely. 

On the continuum level, we therefore have to provide 
equations of motion for each interface point of the crack, 
reflecting either surface diffusion (SD) or the phase trans- 
formation (PT) mechanism between the solid and the gas 
phase. They are coupled to the elastic fields in a nonlo- 
cal and nonlinear manner. The motion is then driven by 
the tendency to lower the total free energy of the system. 
An important and obvious difference is that for SD the 
number of "solid" atoms is conserved, which is not the 
case for the PT mechanism. 

For both transport mechanisms, we consider the 
growth of a single crack in a strip geometry, in order 
to have a constant stress intensity factor. We restrict to 
an effectively two dimensional system by the assumption 
of translational invariance in the z-direction, and assume 
the strip to be infinitely extended in the direction of prop- 
agation, which in our case is chosen to be the x-direction. 
We mainly concentrate on mode I fracture, which means 
that the applied tensile forces act in the y-direction per- 
pendicularly to the crack faces. Apart from this, we will 
also discuss results from the application of mode III load- 
ings, and linear combinations of these two modes. Since 
the crack tip is macroscopically extended, no singularity 
appears and the whole crack pattern can be described 
consistently in the continuum approximation. 

In a Lagrangian description of linear elasticity, the elas- 
tic state of the system is described through a continuous 
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displacement field Ui. Then, the strains are defined as 
the symmetrical spatial derivatives of the displacements, 
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As the total stress field depends linearly on both the 
strain as well as the strain-rate, we conveniently decom- 
pose it into a strain and a strain-rate dependent part, 
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where al*^^ and cr-^"'' are the elastic and viscous stresses, 
respectively, and tik denotes the time derivative of the 
strain Cik. Furthermore, we restrict the considerations to 
fully isotropic media. Then, as given by Hooke's law, the 
elastic stresses are 
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where E is Young's modulus and v the Poisson ratio, and 
we use the Einstein sum convention. By construction, 
the viscous stresses are formally similar to the elastic 
stresses [H], and we therefore write them for a Kelvin 
viscoelasticity model as 
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with the two viscous constants rj and C- 

The evolution of the elastic degrees of freedom within 
the viscoelastic solid is given by Newton's equation of 
motion, and the elastic displacements Ui have to fulfill 
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where p is the mass density. This equation ensures locally 
a force balance between the elastic stress and viscous 
friction on the left hand side and inertia on the right side. 
These equations have to be supplemented by mechanical 
boundary conditions, which are given below for the two 
different transport mechanisms. 



A. Surface diffusion 

For crack growth by surface diffusion, the crack is filled 
with vacuum, and therefore we impose stress free bound- 
ary conditions at the crack contour. 



-I- pVnUi = 0, 



(6) 



where n is the direction normal to the interface, and Vn 
is the normal interface velocity (see Fig. [2]). The second 
term on the left hand side accounts for momentum con- 
servation at the solid- vacuum interface. We point out, 
that in the dynamic limit, when the crack propagation 




FIG. 2: Schematic picture of the steady state crack propa- 
gation by surface diffusion. The crack contour, indicated by 
the solid black line, separates the viscoelastic medium from 
the advancing "vacuum bubble" . During the propagation the 
total amount of solid material is conserved. 



velocity v is of the order of the materials sound speed, 
this term becomes important [3]. 

So far, for an arbitrarily given crack shape and known 
strain history, the mechanical problem is unambiguously 
determined and can be calculated by Eqs. ([l])-(|6| to- 
gether with the outer boundary conditions at the borders 
of the strip, which specify the externally applied loading. 
Next, we have to formulate an evolution equation for 
the crack contour. The motion of the interface is caused 
by thermodynamically induced mass transport processes, 
which diminish the total free energy of the system. The 
local driving force for crack propagation is given by the 
chemical potential p at the solid vacuum interface |42| . 
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with 7 being the surface energy per unit area and k the 
surface curvature, which is counted to be positive, if the 
crack shape is convex; the atomic volume Q appears since 
the chemical potential is defined as free energy per par- 
ticle. We point out that the viscous stresses do not ap- 
pear in the chemical potential, since viscous dissipation 
is a sole property of the bulk, whereas the chemical po- 
tential is needed to describe energy dissipation through 
the motion of the interface. Furthermore, we note that 
due to inertial effects, also the kinetic energy density ap- 
pears in the chemical potential. Counterintuitively, it ap- 
pears with sign opposite to that of the potential energy; 
this can be derived rigorously from variational principles 
[381143]. 

For surface diffusion the motion of the crack surface is 
proportional to the divergence of a flux of solid material 
along the interface. This flux of material is induced by 
gradients of the chemical potential. We express the mo- 
tion of the interface by the local normal velocity Vn and 
obtain 



(8) 



where d/ds denotes the tangential derivative and the dif- 
fusion coefficient Dg has a dimension [Ds] = m^s~^. We 
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FIG. 3: Sketch of a propagating crack, where the phenomenon 
of fracture is interpreted as a phase transformation process 
from a viscoelastic solid to a "dense gas" phase. The crack 
surface, indicated by the solid black line, separates the origi- 
nal viscoelastic medium from the growing "dense gas" phase. 



note that for surface diffusion the amount of solid mate- 
rial is conserved during the crack propagation. A typical 
steady state crack shape using surface diffusion is shown 
in Fig. [2] One can see that the crack first opens up to a 
tip diameter 2/i, and then closes again due to the condi- 
tion of material conservation. 

We point out that this description of mass transport 
is not limited to surface diffusion in its literal sense only. 
Often, many complicated physical processes like plastic 
bulk flow take place in a small zone around the tip. As- 
suming that this zone is relatively thin, the mass trans- 
port can effectively be described by surface diffusion, 
where the detailed information about the process zone 
is hidden in the diffusion coefficient in the spirit of a lu- 
brication approximation. 



B. Phase transformations 

Here we discuss crack propagation by means of a phase 
transformation process, where the solid matrix trans- 
forms into a "broken gas phase" with vanishing elastic 
moduli. We assume that the gas phase and the viscoelas- 
tic medium to have equal mass densities p. Furthermore, 
the interface between this "dense gas" phase and the 
medium is considered to be coherent, i.e. the displace- 
ments are continuous there. With these assumptions, 
two central simplifications are achieved. First, instead of 
Eq. (|6| the mechanical boundary conditions are 



0, 



(9) 



which means that no velocity dependent correction ap- 
pears here, since by the continuity of velocities and den- 
sities also the momentum is contiuous. Second, the ex- 
pression for the chemical potential is replaced by 
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where the kinetic energy contribution does not appear. 
The reason for this simplification is the continuity of the 



kinetic energy density, because the above expression of 
the chemical potential should be more correctly be inter- 
preted as the chemical potential difference between the 
adjacent phases Notice that the inner phase is 

assumed to be infinitely soft, therefore it has a vanishing 
elastic energy density. Again, the motion of the interface 
is locally expressed through the normal velocity, which in 
this case is direct proportional to the chemical potential 
difference at the interface, 



D 

Vn = —fJ. 



(11) 



with 



[D]=m2s-^ 



a kinetic coefficient D having the dimension 
Of course, using this model, the amount 
of solid material is not conserved during crack propa- 
gation. In this sense, our model is strongly related to 
phase field models of fracture based on a non-conserved 
order parameter [44U46j . The crucial difference is that 
the current model is based on well-defined sharp interface 
equations, and therefore the predictions do not depend 
on inherently numerical parameters like a phase field in- 
terface width. A typical steady state crack shape using 
phase transformations is shown in Fig. [3] In contrast to 
surface diffusion the crack keeps its opening, and does 
not close far behind the tip (in the reference state) due 
to the absence of material conservation. 

However, although surface diffusion seems to be more 
adequate for a description of fracture, from a numerical 
point of view the treatment of Eq. ([s]) is much more time- 
consuming due to the higher order spatial derivatives, 
which are not present in Eq. (11). Therefore, modeling 



of fracture as a phase transition process offers numerical 
advantages. 



III. CRACK PROPAGATION: SELECTION 
PRINCIPLES 

The self-consistent selection of the crack velocity, the 
tip radius and even the entire crack shape is a central as- 
pect of the present theory, and will be discussed in more 
detail in this section. The bulk equation ([s]), in combi- 
nation with equations (|6|-([8]) for the surface diffusion, or 
equations ([9])-(ll| for the phase transformation mecha- 
nism, describes the dynamics of the two models. In both 
cases they lead to a complicated free boundary problem. 

Before starting to solve the full free boundary problem 
numerically, we first discuss qualitatively the existence of 
steady state solutions, by the use of scaling arguments. 
Here, the term steady state describes a non-equilibrium 
solution, for which the crack is moving with a constant 
velocity v, and in a co-moving frame of reference - fol- 
lowing the crack tip with the same steady state velocity 
V - the shape is constant in time. The following scaling 
arguments are fairly generic and can similarly be applied 
to both the SD and the PT model and predict the char- 
acteristic velocity and tip scale. 
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We will address the selection problem on different lev- 
els: First, we use pure dimensional arguments for po- 
tential length- and velocity scales, and we show that 
only with the additional parameters stemming from vis- 
coelasticity or mass density microscopic tip scales appear, 
which can select a tip radius. This argument captures al- 
ready the essential physical situation, and therefore the 
following, more extended discussion could be skipped for 
first reading. There, we revisit the behavior with a more 
detailed analysis of the equations of motion; the basic 
outcome of this more advanced investigation is that for 
the purely elastostatic case tip radius and velocity cannot 
be selected inpendently. We then show how the inclusion 
of inertial or viscoelastic effects cures this problem. 

The simplest assessment is to determine the possible 
nontrivial parameter combinations in the model in order 
to form a lengthscale. For pure linear static elasticity, 
however, it is not possible to set a microscopic length 
scale out of the material parameters and the applied load, 
and therefore selection is not possible in this case. The 
lack of such a lengthscale is the reason for the cusp sin- 
gularity of the ATG instability and the impossibility of 
a steady state crack growth under these conditions. 

The situation is different, when inertial effects are 
taken into account, because then the Rayleigh wave speed 
(the sounds wave speed at a free surface), vr ~ {E/pY^^ 
appears as characteristic velocity scale in the problem. 
Then we expect the crack velocity to be of this order, 
vr, and we can form a microscopic lengthscale by 
the parameter combinations {Dg/vnY fo^' ^-D ^-^d D jv^ 
for FT, which we expect to be the tip scales in this 
case. If instead of inertial effects viscous bulk damp- 
ing becomes relevant, the characteristic velocity scale 
is Vq - {D.E^ jrf'fl^ for SD and ~ [DE/rjf/'^ for 
FT. Consequently, we expect the characteristic tip scales 
{Dsv/Ey/^ for SD and {Dr)/Ef/'^. 

For a more detailed analysis, we have to inspect the 
fields and the equations of motion. We first return to 
the case without viscous or inertial effects. We note that 
the stresses on the boundary of the crack tip with finite 
radius tq scale as 



Kr, 



-1/2 



(12) 



where K is the stress intensity factor. Since this anal- 
ysis can be applied both to mode I and mode III, we 
do not indicate the loading conditions in the stress in- 
tensity factor. Also, by boundary conditions the normal 
and shear stresses on the crack surfaces vanish (or are 
small, if the momentum transfer term is included for sur- 
face diffusion, provided that the crack speed is substan- 
tially smaller than the speed of sound). Therefore, the 
only nontrivial stress component is the tangential stress, 
which depends on the shape of the entire crack. 

Although the stress scaling ( 12 1 appears to be natural 



in the framework of fracture mechanics, it is not trivial, 
and therefore we discuss it in more detail. Since we intend 
to describe cracks with a finite tip radius, the stress field 
typically contains not only r~'^/'^ terms but also faster 



diverging terms 7--3/2^ ,,-5/2^ ^ _ _ These terms cannot be 
present for sharp crack tips, since they would lead to a 
diverging elastic energy, but cannot be excluded for finite 
ro, since then the divergency is cut off. Therefore, the 
stress field typically consists of singular and regular parts 
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K 
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(13) 

in a polar representation. For the sake of brevity we 
do not write the angular dependence /^*' explicitly. The 
regular part of the stress Greg contains only constant con- 
tribution ("T stress") and positive powers of r. In the 
region rp < r <C W, i.e. close to the (small) crack tip on 
the scale of the system size W , the ascending powers can 
however be ignored. Also, the T stress scales as K/W^/'^, 
and therefore vanishes for large system sizes and given 
stress intensity factor K. Consequently, areg = for 
the further discussion. The above representation is the 
heart for the multipole expansion method that will be 
introduced below in Section |IV A[ Notice that the higher 
order modes seem to violate the anticipated r^J"^ scaling 



of Eq. (12). To understand this situation, we inspect a 



crack with a different tip radius, which we obtain by a ge- 
ometrical rescaling r — > r/a, and with the same stress in- 
tensity factor. We make the scaling ansatz for the stresses 
a for this rescaled geometry as a{r,9) = f3a{r/a,9) 
in relation to the original problem. At large distances 
from the tip, the main mode a ~ if/(27rr)^/^ prevails, 
hence a{r,e) = Pa{r/a,e) ~ Pa^/^K/{2TTry^^. On the 
other hand, both cracks look identical on large distances, 
where they become sharp straight cracks, and therefore 



a 



-1/2 



On the crack surface, where also the higher 



modes are relevant, the stress of the rescaled problem is 
therefore a{fo,9) = (7{aro:0) = a(ro,9)a~^/'^ . Hence, 
for a crack with a four times sharper tip, the surface 
stresses are two times higher, as stated in Eq. (12). 



Obviously, the interface curvature in the tip region 
scales as K ~ l/j'o- Hence, as long as only static lin- 
ear elasticity is taken into account, all contributions to 
the chemical potential scale as /x ~ ■ Consequently, a 
rescaling of the equations of motion ([s]) or (11) is pos- 
sible: Formally, the equations of motion depend only 
on the dimensionless combinations vr^/Ds for the SD 
mechanism and vr^/D for the FT dynamics. All other 
parameters combine to the dimensionless driving force 
A — K^{1 — v'^)/2E^ (in case of pure mode I loading), 
where A = 1 corresponds to the Griffith point. In other 
words, the radius rg and the steady state velocity v can- 
not be selected separately within the framework of the 
pure static theory of elasticity. Even if a steady state so- 
lution exists, it would still degenerate to a one-parameter 
family of solutions with either fast cracks with small tip 
radii or slow cracks with large tip scales. It turns out, 
however, that no steady state solutions exist, which is 
exactly the aforementioned cusp singularity of the ATG 
instability. 

Additionally to this inspection of the behavior in the 
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tip region, we can get more insights from the analysis 
of the crack shapes in the tail region, where the elastic 
stresses have decayed. To that end we assume that the 
crack is growing in the steady state regime in positive x- 
direction with a constant velocity v. For the SD model, 
the shape equation ^ can be integrated once, and in the 
co-moving frame of reference we obtain 

7^2 ds 

This is a complicated, non-linear third order equation 
with non-local contributions arising from the elastic 
fields, since aik depends on the entire shape. In the 
tail region the shape equation is simplified to the third 
order differential equation Dgy'" — vy due to the de- 
cay of the stress fields, which has two growing and one 
decaying solution. Only the latter y{x — >■ — cxi) = 
Ae'xjp({v/ Ds^^^^x) asymptotically describes the physi- 
cally allowed shape. We switch to a polar coordinate 
system x = r{9) cos9,y = r{9)sa\9 for the tip region, 
and focus on symmetrical solutions, r{9) = r{~9). Since 
the physical properties, curvature and stresses, do not de- 
pend on the choice of the coordinate system but only on 
the crack shape, we can arbitrarily chose r(9 = 0) = tq, 
with the a priori unknown tip radius = 1/k(0). Then 
from symmetry considerations and the definition of the 
tip curvature, k = (r^ -|-2r'^ — rr")/(r^ -t-r'^)^/^, the nat- 
ural conditions ^'(O) = r"(0) = arise. Integration over 
the upper interface > requires the suppression of the 
two growing exponentials at the tail, which imposes two 
boundary conditions. For a given external loading, these 
two conditions can be fulfilled by a proper selection of 
the tip radius ro and growth velocity v. However, since 
the use of only linear static theory of elasticity does not 
allow the independent selection of both the tip radius tq 
and the steady state velocity v, the selection will not sup- 
press both growing exponentials at the same time, and 
consequently a crack like solution does not exist '35]. 

For the FT model, a similar argument can be given 
[36t l37] . In the tail region, where the elastic stresses 
have decayed, the shape equation becomes simply —vy' = 
Dy" . Its general solution, y{x) = h + B eyi\i{—vx/ D), 
contains the finite crack opening h and a growing expo- 
nential. Notice, that in contrast to the surface diffusion 
process a finite opening 2h cannot be excluded since we 
do not have to obey mass conservation here. Suppress- 
ing the exponential and selecting a finite tail opening h 
finally requires again the independent selection of both 
the steady state propagation velocity and the crack tip 
radius. Consequently, a steady state solution for a grow- 
ing crack in the framework of the static theory of elastic- 
ity does not exist. 

The situation changes if additional length scales en- 
ter into the description, and two natural aforementioned 
extensions for a description of crack growth are viscous 
bulk dissipation and dynamic elasticity. Nonlinear elas- 
tic or plastic effects are also conceivable, but are beyond 
the scope of the present article [331 . 



Although the viscous stress defined in Eq. Q intro- 
duces two new parameters, the timescales, which are in- 
troduced by them, should typically be of the same order. 
By setting = v we restrict to the case of only one ad- 
ditional time scale r — rj/E to simplify the situation. 
Then, considering static elasticity and viscous bulk dis- 
sipation, additionally the dimensionless ratio v/vq with 
vq = (plTfl'^ for FT and vq = (D.T-^y/^ for SD ap- 
pears in the equations of motion, and therefore a rescal- 
ing is no longer possible. Then this additional free pa- 
rameter allows to independently select both the tip scale 
ro and the steady state velocity v, so that the two grow- 
ing exponentials can be suppressed. 

To make this more explicit, we note that by virtue 
of Eqs. (jsl) and ^ for steady state growth cr^^'^^ = 
—vrdxCr'^^ . Therefore, as consequence of the force bal- 
ance condition ([5| we get a correction to the elastic 
stresses which depends on the dimensionless parame- 
ter vt/vq. For SD, in the interface evolution equation 
the two non-dimensional parameters Si = vr^/Dg and 
S2 = VT Itq appear, which contain different combinations 
of the tip radius and the crack velocity, hence these two 
nonlinear eigenvalues can be selected independently to 
suppress the two growing exponential terms in the tail re- 
gion. Since s\ and S2 are of order unity for driving forces 
of order one, we therefore get s\s\ = v'^T^/Dg ~ 1, hence 
V ^ {DsT~'^y^^ , which is the predicted velocity scale vq. 
Analogous arguments can be used for the FT mechanism. 
Similarly, we obtain v ^ vr if inertial effects are relevant, 
and the microscopic tip scales are rg ~ {Ds/vnY for SD 
and ro - D/vr for FT [Ml EH HH] . 

Altogether, we conclude that independent of the con- 
sidered mass transport mechanism, steady state growth 
of cracks is possible if apart from static elasticity at least 
one additional effect is taken into account. Furthermore, 
for both mechanisms a tip splitting is at least conceiv- 
able for high applied tensions due to a secondary ATG 

— 1/2 

instability: Since a ^ Ktq in the tip region and the 
local ATG length is Lq ^ E'j/a'^, an instability may oc- 
cur, provided that the tip radius reaches by blunting the 
order of the ATG length. 

The similarity of the scaling arguments to predict 
steady state growth for both SD and FT emphasizes 
the close connection between the two models. From a 
physical point of view the SD model is probably more 
appealing, but also more difficult to solve numerically. 
However, the preceding arguments suggests that many 
generic properties of the model should also be reflected 
in the simpler FT model. In section [V] we will give a 
detailed comparison between the model behaviors. 



IV. NUMERICAL METHODS 

The free boundary problem, which arises from the cou- 
pling of non-local dynamic elasticity or viscoelasticity to 
the interface kinetics, is studied by the use of two com- 
plementary methods, which are presented in this section. 



8 



The first method is a sharp interface method, based on 
the expansion of the elastic fields in a series of eigenfunc- 
tions of a straight mathematical cut. This multipole ex- 
pansion method, designed to simulate efficiently steady 
state crack propagation, delivers precise results in two 
limiting cases: Slow crack propagation with viscoelastic 
effects and elastodynamic fracture without bulk dissipa- 
tion. Situations where both effects play a role can only 
be treated in a perturbative manner. 

As second method we use a fully dynamic phase field 
model with a sharp interface limit corresponding to our 
model equations. In contrast to the multipole expansion 
method the phase field approach allows also to investigate 
transient behaviors and crack branching, and also enables 
to model both elastodynamic and viscoelastic effects in 
a uniform framework. However, obtaining quantitative 
results comparable with those by the multipole expansion 
method is computationally very expensive. 



A. Multipole Expansion Method 



finding proper expansion coefficients Ci in order to satisfy 
the boundary conditions Eqs. ([6]) or ([9]). This reduction 
makes this method numerically very efficient. 

However, to our best knowledge, there is no closed so- 
lution to the full problem, which means that the eigen- 
functions for the visco-elastodynamic problem of a mov- 
ing mathematical cut are not known. Therefore, we focus 
here on two limiting cases of Eq. (14 1: First, the static 



limit of viscoelasticity, where v vji, and therefore the 
term on the right hand side is neglected, and second the 
elastodynamic limit where the viscous damping vanishes, 
i.e. r = 0. Here, we mainly deal with mode I fracture, 
and therefore we illustrate the corresponding procedures 
to solve both the viscoelastic problem and the elastody- 
namic problem for these loading conditions. The subse- 
quent technical subsections for the determination of the 
eigenmodes for these two cases can be skipped for the first 
reading. For mode HI loading, which is mathematically 
simpler, similar approaches can be found; in particular, 
the solution of the viscoelastic mode HI problem is pre- 
sented in l39l. 



For solution of the steady state problem of crack prop- 
agation with the multipole expansion method we divide 
the problem into two parts: First, the solution of the 
mechanical problem for an arbitrary, but known crack 
shape and velocity f, and second the evolution of the 
crack contour and adjustment of the velocity. These two 
steps are iterated until a self-consistent steady state so- 
lution is found. 

To simplify the appearance of viscosity in our equa- 
tions, we assume = v and thereby focus on the case 
of only one additional time scale t = t]/E due to viscos- 
ity. With this simplification the dissipative stress tensor 
is related to the elastic stress tensor by a. 



{vis) 
ik 



.lei) 
'ik 



We note that for mode HI fracture such a simplifying pa- 
rameter choice is not necessary, since there always only 
one timescale appears. For steady state growth, the time 
derivative in the co-moving frame is then replaced by a 
spacial derivative with respect to the crack propagation 
direction x, d/dt = —vd/dx. Consequently, the steady 
state mechanical bulk equilibrium equations Eq. ([5| , con- 
taining both viscoelastic and inertial effects, become 



d 
dxk 



Ml) 



d (el)\ 2^""" 



dx"^ 



(14) 



The basic idea for solving the elastic problem, is to 
write the elastic fields as an expansion in eigenfunctions 
of the differential operator corresponding to the equa- 



tions of motion (14) for a straight moving cut. Formally, 



the structure of the stress fields becomes 



K 
(27rr)i/2 



1 + Ci 



/(I) /(2) 



C2- 



(15) 



where the coefficients Ci are the amplitudes of the eigen- 
modes f^'^\ Then, the bulk equations are automatically 
fulfilled, and the problem is reduced to a linear one for 



1. Viscoelasticity 

First, we consider the limit of small crack velocities, 
i.e. V <^ vr. In this limit the the term from inertia on 



the right hand side Eq. ( 14 1 can be omitted. Therefore, 



the force equilibrium condition in the static limit for the 
steady state situation reads: 



d 

dxk 



Ael) 
'ik 



d (el] 



= 0. 



(16) 



For the solution of this problem, we use Airy functions, 
which are well known in static elasticity. Here, we gen- 
eralize this approach to viscoelastic materials. 

We begin with the treatment of a static elastic problem 
and introduce for convenience a complex Airy function 
U{z), with z — x + iy. The usual (real) Airy function is 
defined as its real part. 



U{x,y)^n{U{z)). 
The usual relations to obtain stresses are 



dxdy ' 



dx-^ 



(17) 



(18) 



Compatibility, i.e. the existence of a displacement field 
from which the elastic strains can be derived, is equiva- 
lent to 



A^U = 0. 



(19) 



In most cases, the complex Airy function U is not ana- 
lytic, and the reason is that its real part has to satisfy 
only the biharmonic equation and not the Laplace equa- 
tion. We therefore make the ansatz 



U^f + zg, 



(20) 
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with f{z) and g{z) being analytic functions (apart e.g. 
from a branch cut for crack problems); the bar denotes 
complex conjugation. This means that with / = /i +1/2 
and real functions fi{x,y), f2{x,y) the Cauchy-Riemann 
equations hold: 



dx 



dh 
dy ' 



dh 
dy 



dh 
dx 



(21) 



With the above structure Eq. (|20|) the biharmonic equa- 



tion Eq. (191 is automatically fulfilled. Stresses can be 



expressed as 



<y.. = ^[-r + 2g'~zg" 
c^xy ^ +zg \, 



(22) 
(23) 
(24) 



where the ' denotes the derivative with respect to z. The 
expression for displacements are 



2/i(u:r + iuy) = (3 - 4i^)g - zg'{z) - f'{z), 



(25) 



with /i = E /2{1 + v), and thus we get for the derivatives 
and strain componentents 



dux 
dy 

dUy 

dx 

f-xy 



1 

2^1 

1 

2^1 
1 

241 
1 

241 
1 

241 



3?[2(l-2z.)g'-z5"-/"], 
3?[2(l-2i.)g' + z5" + /"], 
^[-A{l-v)g' + zg" + r'], 

^m-v)g' +zg" +n, 



3 [zg" + /" 



(26) 
(27) 
(28) 
(29) 
(30) 



As mentioned before, the total stress decomposes into 
an elastic and a viscoelastic contribution, 



cr. 



(tot) 



Ml) 



ik ^ "ik ^ik ' 

where the latter is for steady state growth 



{vis) 



(vis) 



d (el) 



(31) 



(32) 



and consequently we have the force balance condition 



(16). In principle, it is not required that elastic and 



viscous stress satisfy the force balance separately, but 
only Eq. (16 1 must hold for the total stress. Also, only 



the elastic fields needs to satisfy compatibility conditions. 
However, as we will see, all fields fulfill these conditions 
even separately. 

The ansatz is that both the elastic and the total stress 
field can be derived from Airy functions which satisfy 
the biharmonic equation. In particular, we anticipate 
the following structure of the complex Airy functions 



(33) 



with analytic functions f{z),g{z),F{z),G{z); for the 
present crack problem these functions are not analytic 
everywhere, but have a branch cut along the negative 
real axis, see below. As we have seen above, the real 
Airy functions [/(*°*) = 3?(Z^(*°*)), U'^"^^ = 3?(Z^(^')) then 
satisfy the biharmonic equation. Stresses can be derived 
from Eqs. (p2l)-((24l). 



Provided that the following equation is fulfilled, 



VT—U 

dx 



(el) ^ jjitot) 



(34) 



then the steady state force balance ( 16 1 is fulfilled and 



a valid elastic displacement field exists by construction. 
We note that the complex Airy functions are not diffcr- 
entiable in the complex sense due to the appearance of 
the complex conjugate factor z, and thus we cannot gen- 
eralize dx^ — ^d/dz. However, with the above ansatz 
([33| the equation (pll) is fulfilled if 



f + zg-vT{f + g + zg') = F + zG (35) 
holds. Separating "harmonic" and "biharmonic" parts 



gives 



f-vrif + g) = F, 
g - vTg' = G. 



(36) 
(37) 



Again, if (36) and (37) are satisfied, then (35 1 is also 
valid. 

We write the functions F and G now as series expan- 
sions in the set of eigenfunctions of a straight mode I 
cut. 



m— — 1 

00 

G = ^ Br^z^l^-'-. 

?n— 



(38) 
(39) 



Notice that the summations start from different values 
of TO, because the function G appears with an additional 
prefactor z in the complex Airy function. The lowest 
value of m corresponds to the main mode. In order to 
have the correct mode I symmetry, the coefficients of ex- 
pansion Am and i?m are real. 

The far field behavior is controlled by the term with 
the lowest value of to, which is the main mode. On large 
distances r from the tip, the crack looks like a semi- 
infinite mathematical cut, and this is refiected by the 
proper r^^/^ decay of the stresses in this purely elastic 
regime. The prefactor of the main mode is therefore re- 
lated to the stress intensity factor, thus 



A. 



Kr 



3V2n 



(40) 



Also, wc have the requirement that on the straight cut 
normal and shear stresses have to vanish, hence 



Bo = 3A^i. 



(41) 
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We write now also the functions / and g as series 

oo 

f ^ J2 (42) 

BrnQm (43) 



OO 

E 

m=0 



with analytical functions fm and gm, which will be de- 
termined below. We define 



4 — 



Bn 



(Am ^ 0), 
{A,n = 0). 



Notice that we assigned for convenience 
B_i = 0. 

Provided that 

l/2-m 



for m = 0, 1, 2, 

fin 

for A„i 7^ or 



5m - vrg^ = 
., and either 

A 



f 



(45) 



(46) 



(47) 



(48) 



for Am — hold for m ~ —1, 0, 1, . . ., then Eqs. (36) and 
(371 are satisfied. 



Notice that the distinction between the regular case 
Am 7^ and the singular case Am = is relevant also for 
practical purposes: For numerics we cut off the expansion 
of F at M - 1 and G at M ; then Bm 0, but Am = 0, 
so for the last mode we always encounter this situation. 



Obviously, the equation for / (47) can be solved as 



soon as the solution of the equation for g ( 46 ) is known 
;ene 



The homogeneous equation for g, g^ —vrg^m = has 



the solution gm' = Dm exp(2:/?;T). Variation of constant 
Dm — >■ Dm{z) then gives the general solution of Eq. (46) 



gm = Dmiz) exp{z/vT) 



(49) 



with 

Dm{z) 



1 

VT 



exp 



{~z/vT)z^''^~"'dz + const, (50) 



where the integration constant has to be chosen such that 
exponential growth terms in ( 49 1 are suppressed. We ob- 



tain in particular for m = — 1 with the proper integration 
constant 



I?_i(z) — exp(— z/wr) 



^3/2 



+ 7:VT z 



1/2 



(i;t)'^/^ V^erfci/z/wr 



(51) 



with the (complex) complementary error function erfc. 
Thus we obtain 



9-i{z) 



^3/2 



-|- -UT Z 



1/2 



+ -{vt)'^/'^^/tt eyiY>{z/vT)eT:ic^J z/vT. (52) 

All higher modes can be obtained from the recursion re- 
lation 



1 



9m+l 



(ij — m)vT 



^1/2-m 



(53) 



(44) which follows from Eqs. (|49| and (50) and the proper 



choice of the integration constant. In particular, 

50 = z'^f'^ + ^(uT)^/^v^exp(z/ur)erfcA/^:/wT, (54) 
which is the main mode term. 



The equation for /, Eq. (47), is treated in the same 
way, and we obtain 



fni — 



3m + -T^^T" hm {Am ^ 0) 
[Am = 0), 



A. 



(55) 



with another analytical function hm{z). It obeys the re- 
cursion relation 



1 



and 



m-t-l 



- ,3/2 



Qm ~1~ ^m] 



(56) 



15 

T 

5 



.1/2 _ 



VT z"'' r(wT)^^^\/7i' X 



VT exp 



{z/vT)eTicy'^^JvT. (57) 



In particular 
ho 



/2 



2 

X I z 



-(vt) 



-1/2. 



VT exp 



[z I VT^QTic^fzJx 



(58) 



To summarize, Eqs. (53)-(57) provide a series of 



eigenfunctions for the steady state equation of motion 
Eq. (16). From these eigenfunctions, via Eqs. (42) and 



(43), together with (22)-(24) the total stress field can be 



calculated as a function of the coefficients of expansion 
Aq, Ai, . . . and Bi, B2, ■ ■ ■■ While the main mode coeffi- 



cients A_i and Bq are given by Eqs. (40) and (41), the 



remaining coefficients {Ai}, {Bi} are now determined in 
order to fulfill the conditions cr„„ = cr„s = on the actual 
crack contour (n and s are normal and tangential direc- 
tions respectively) . The optimization of these expansion 
coefficients is equivalent to finding the minimum of the 
function 



Ri{A,},{B.}) 



2 

nn 



ds 



(59) 



with respect to {^i}, {Bi}., where the integration is per- 
formed along the crack contour. 
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2. Elastodynamics 

Now we discuss the solution of the elastic boundary 
value problem in the dynamic limit of vanishing viscous 
bulk dissipation, i.e. r = 0. Therefore we briefly review 
the analysis given in . Following Ref . [H [5T] , we in- 
troduce two real functions (j){x,y^t) and ip{x,y,t) which 
are related to the displacements as follows, 



where f[k\{^i ^) ^''^d fl']^'a{0, v) are the universal angular 
distributions for the dilatational and shear contributions 
which also depend on the propagation velocity. In anal- 
ogy to the procedure above, the unknown coefhcients of 
the series expansion are determined by minimization of 
the residuum 



d(j) dip 
dy 



dx dy ' 

Using the decompositions of the displacement field, the 
steady state bulk equations (14 1 become homogeneous 
Laplace equations 



dx"^ dy\ 



where the coordinates perpendicular to the crack are 
rescaled by either yd — aay or j/s — Ugy for either 
the function (f> or the function ip. Here, we have de- 

= 1 — ij'^/c^, where, Cd ~ 



fined = 1 — jc^ and a 
^E{\-v)Ip{\-2v){\^v) and = sjEl2p(\ -h v) 
are the dilatational and shear sound speeds respectively. 
Since we are looking for solutions obeying the mode I 
symmetry, we propose the ansatz 



E/i 3/2-n 
A,^r^ cos 2 

n=0 



00 



B r^l'^~ 



(61) 



(62) 



in rescaled polar coordinates, which are related to the 
co-moving Cartesian coordinates via a; = cos 6d = 
Ts cos 9s, yd = Td sin 6d and j/s — sin Og . For a crack 
with a sharp tip, only the mode with n = is allowed, 
which corresponds to the usual a ^ r^^/"^ behavior. For 
this mode, the boundary conditions on the straight cut 
and the matching to the far field behavior demand 



Bo 



where Kf^" is the dynamical mode I stress intensity fac- 
tor [3], related to the static stress intensity factor as 

.2^2 \ 1/2 



2TT3E{4:asad 
2ad 



a 



2\2\ i 



dyn 



1 



:Ao 



at 



(63) 
(64) 



K 



dyn 



j^stat 



(1 



^ 4:asad - (1 + a ^sY 
ad{l - Q;2) 



(65) 



Each eigenmode of Eqs. (61) and (|62j) satisfies the elas- 
todynamic bulk equation (60). The coefficients Aq and 
Bo are determined by Eqs. (63) and (64) for the correct 



far-field behavior, whereas all other modes decay faster. 
Consequently, we obtain the formal stress field expansion. 



dyn 



(27rr)i/2 



/ N=oo A f(") I R f(»*) 

^(0) ^nJik,d ^ ^nJik^s 

J ik ' / , 

K n=l 



(66) 



i?({A.},{Bj) 



{cTni+PVnUi) ds 



(67) 



with respect to the coefficients A^ and Bi, for a given 
crack contour and steady state velocity; the integration 
domain is the crack contour. Notice that this residuum is 



used for SD, whereas for FT it is the same as in Eq. (59 ) 



3. Approximative viscoelastodynamic model 

To our best knowledge, there is no exact solution of 
the full problem given by Eq. (14), which contains both 
dynamical and viscous effects. We therefore suggest an 
approximative model, which captures essential physical 
aspects and gives exact results both in the limit of van- 
ishing viscous damping and treats viscous damping in the 
sense of a rigorous perturbation theory for low velocities. 

To motivate the model let us first consider the case of 
static elasticity, where inertial effects can be neglected. 
Thus, we solve the elastic problem consisting of three in- 
gredients, i.e. bulk equilibrium da'^f'^^/dxj — 0, stress 
free boundaries at the crack surfaces, a^'l'"'^^ = and 



Kjr at large distances, and additionally 



^(ei,0) 

compatibility, which means that the strain tensor, which 
is related to the stress via Hooke's law, can be derived 
from a displacement field. We use here the superscript 
(e/, 0) to indicate that we are dealing here with a purely 
elastic field, which is later used as zeroth order for a per- 
turbative treatment. By these requirements the solution 
for u^j'"'^^ is formally uniquely defined (apart from trans- 
lational and rotational degrees of freedom). 

On the other hand, for the viscoelastic steady state 
problem the total stress consists additively of the elas- 

'^T^\ where the 



tic and viscous stress, crf'°*^ = alf^ 
viscous stress is related to the elastic stress by cr. 



(vis) 



—vrdxcrlj (we assume here again that we have only one 
viscoelastic timescale and that the crack moves in posi- 
tive X direction). Now the total stress has to fulfill me- 
chanical equilibrium, daf°^^ /dxj = 0, the total normal 
and shear stresses have to vanish on the crack surfaces, 
fjf"^^ = and far away (where the behavior is anyway 
purely elastic) we have the same asymptotic behavior 
ct[*°*'' ^ Kir^^/^ . Also, as we have seen in the section on 
the multipole expansion method, also the totoZ stress field 
satisfies compatibility. Therefore, by uniqueness of the 
solution, the solution for the total stress is here exactly 
the same as before for the purely elastostatic problem. 



(tot) 



(el,0) 
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Now we use the solution of the elastic problem and 
introduce viscosity perturbatively, where we use r as 
small expansion parameter. The zeroth order solution 
'^if'^^ generates a viscous stress to first order in r, i.e. 



(vis,l) 
J- ■ 



(ei,0) 



However, up to first order order 



the sum of these two terms, c-^''^' + cr-J'"'^^ does not yet 
satisfy boundary conditions on the crack surfaces (but 
they do satisfy the bulk force balance conditions), and 
therefore another elastic correction term must appear to 

first order, c^^^^'^K Hence, up to first oder the total stress 

(e/,0) 



IS a 



{tot} 



fJ, 



(ei,l) 



{vis A) 
T- 



-0{t^). On the Other 
hand, by the aforementioned uniqueness of the solution. 



{tot) 



clj"' = (^ij ' ' to all orders. Consequently, we obtain 

{el,l) {vis.l) , o (ei,0) 



■ij " I] ' " ■ "-^" ij 

For the equations of motion we need an expression for 
the chemical potential, which depends on the elastic part 
of the stress only. Hence we get a first order correction 
to the chemical potential 



(1) 



(o) H- 
dx 



(68) 



4- Steady state crack growth 

Once we can solve the mechanical problem for arbi- 
trary shape, we can solve the free boundary problem for 
the steady state crack propagation. The latter is de- 
scribed, depending on the mechanism of propagation, by 
the set of Eqs. in case of the PT model or by 

Eqs. ([7])-([8| in case of SD. The strategy for solving the 
problem is as follows: For a given guessed initial crack 
shape and velocity, we determine the unknown coeffi- 
cients An and i?„ from the boundary conditions. After- 
wards, we calculate the chemical potential and the nor- 
mal velocity at each point of the interface. The new shape 
is obtained by advancing the crack according to the lo- 
cal interface velocities. This procedure is repeated until 
the steady state is reached, which means that the shape 
of the crack in the co-moving frame of reference remains 
unchanged [55]. This "quasi-dynamical" approach pro- 
vides a natural way to solve the problem, as it follows the 
physical configurations to reach the steady state. Then 
the following relation between the local normal velocity 
and the steady state velocity v holds 



0, 



(70) 



Here we have made use of the property cr,|^''°^e|j'' — 

^^ik '^'^ ^fk ^ ■which follows from Hooke's law. Notice that 
there is no need to decorate the strain with a superscript 
(eZ), since strain is by definition an elastic property (the 
viscous stress is related to the strain rate). This descrip- 
tion is a rigorous perturbative treatment of the viscoelas- 
tic theory in the quasistatic limit, v <^ vr. 

However, this concept cannot strictly be extended to 
the case of dynamical elasticity, since there on the right 
hand side of the Newtonian equation the acceleration 
term pili contains only elastic displacements (so even in 
the force balance equation for the total stress the right 
hand side contains the elastic accelerations), and there- 
fore the purely elastodynamic and the total viscoelasto- 
dynamic stresses do not obey the same equations. 



Instead, we use the above recipe (68) to incorporate 



viscous damping in the chemical potential also with in- 
ertial effects, and consider this as an approximative vis- 
coelastodynamic model; of course, this model is not rig- 
orously derived, but captures essential aspects of the 
physics and is still exact for r = and becomes a rigorous 
perturbation theory for v <C vr. Thus, in the framework 
of this model the chemical potential for the PT model is 



where is the x component of the normal vector point- 
ing into the solid phase. This is a purely geometrical re- 
lation and therefore it is independent on the mechanism 
of crack propagation, i.e. it is valid for the SD model as 
well as for the PT model. This equation gives us an al- 
ternative approach to the "quasi-dynamical approach" . 



Namely, we directly solve the nonlinear equation ( 70 ) as 



a functional of the crack shape and the tip velocity v 
by Newton's method complemented by Powell's hybrid 
method [53, ,54^ and we refer to this as the "steady state 
approach". We stress here that the "steady state ap- 
proach" is preferable especially in case of the SD model, 
where we thus can avoid to solve the fourth order differ- 
ential equation ([s]). 

Finally, we define the dimensionless driving force 



A = A/ + A 



/// 



-Kj H 



2Ej 



2Ej 



iii^ 



(71) 



where we also include the possibility of mixed- mode load- 
ing. Here, A = 1 corresponds to the Griffith point, and 
the energetics necessarily require A > 1 for crack growth. 



B. Phase Field modeling of fracture 



fj, = fl 



A^{0) (0) 
2 ik ik 



,(0)^ 

''^ dx 



7K 



(69) 



where the stresses a'^^} and strains e^-^'' are calculated 



from the elastodynamic eigenfunctions (|6l[)"(62 ). No- 



tice, that for the SD mechanism, we also have to account 
for the kinetic energy density as in Eq. 0. 



During the past years, phase field modeling has 
emerged as a promising approach to model crack prop- 
agation by continuum methods (see |26| for a recent re- 
view). This method is especially advantageous due to its 
high versality to study quite complicated crack patterns 
as well as multi crack situations [55]. Nowadays, phase 
field models capture many known features of cracks |441 - 
1311 [SS] ; However, a significant attribute of most of these 
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descriptions is that the scale of the growing patterns is 
always set by the phase field interface width, which is a 
purely numerical parameter and not directly connected to 
physical properties; therefore these models do not possess 
a valid sharp interface limit. Alternative descriptions, 
which are intended to investigate the influence of elastic 
stresses on the morphological deformation of surfaces due 
to phase transition processes, are based on macroscopic 
equations of motion. But they suffer from inherent finite 
time singularities which do not allow steady state crack 
growth unless the tip radius is again limited by the phase 
field interface width [551150] . 

Since the phase field method was originally developed 
to mainly simulate the dynamics of solidification pro- 
cesses, it is of course more natural to formulate a phase 
field model for fracture using the PT mechanism Eq. (111. 



However, we mention here that it is also possible to for- 
mulate phase field models for crack propagation by sur- 
face diffusion [571 - [5^ . and for example the initial stage 
of the ATG-instability has already been reproduced us- 
ing such kind of phase field models. Nevertheless, for 
the current purpose, we restrict the discussion to phase 
field modeling for crack propagation using non-conserved 
order parameter dynamics, see Eq. (11). 



For the formulation of the present phase field model, 
we start with the introduction of a continuous phase field 
(/), which will discriminate between the different material 
states. We define </) = 1 for the viscoelastic medium, 
and (p — ior the "broken phase". This region does 
not support clastic stresses (the material is broken), but 
still it has the same density as the surrounding matrix. 
Therefore, we use the notation of a "dense gas phase" to 
underline that we do not have vacuum inside the crack. 

We start from a free energy functional, similar to [50] 



ifs+fd^+fel)dV, 



(72) 



where fsC^c/)) — 37^(V0)'^/2 is the gradient energy den- 
sity and fdw{4>) = 670^(1 — 4>Y /S, is the double well po- 
tential, guaranteeing that the free energy functional has 
two local minima at = and cj) — 1 corresponding to 
the two distinct phases of the system. The form of the 
double well potential and the gradient energy density are 
chosen such that the phase field parameter ^ defines the 
interface width and the parameter 7 corresponds to the 
interface energy of the sharp interface description |58| . 
Finally the elastic energy density contribution is 



fe 



h{(j))E 

2(1 + !/) 
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(73) 



where h{(j)) — (fP{3 — 20) interpolates the elastic modu- 
lus between zero for the dense gas phase and one for the 
viscoelastic medium. It is the simplest polynomial sat- 
isfying the necessary interpolation conditions h{0) = 
and h{l) = 1 and having a vanishing slope at = and 
= 1, in order not to shift the bulk states. Here, for the 



sake of brevity we consider the Poisson ratio to be phase 
independent. 

The evolution of the elastic fields is determined by the 
principle of momentum conservation according toEq. ([5|, 
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(74) 



where the elastic stresses are defined as the derivative of 
the elastic free energy density with respect to the strains. 



I.e. cr. 



(el) 
ik 



dfei/dtik- In order to have vanishing viscous 
damping inside the crack, we use define the viscosity pa- 
rameter 77 to be phase field dependent, i.e. 77 — rih(4>) in 
Eq. (|4]), while the parameter C remains phase indepen- 
dent. 

The phase field dynamics is related to the functional 
derivative of the free energy with respect to the phase 
field variable, 



90 _ D rsF 



(75) 



Ui— const. 



where D corresponds to the above mentioned kinetic co- 
efficient of the phase transformation model. Here it is 
assumed that the viscous dissipation does not affect the 
phase field dynamics. According to our sharp interface 
model of crack propagation, we consider viscosity to be a 
bulk property, which does not affect the phase change be- 
havior directly. We ignore local heating effects through 
bulk or interfacial dissipation, assuming that the heat 
diffusion is sufRciently fast. 

Using the phase field method, we investigate crack 
growth in a strip geometry with fixed displacements at 
the upper and lower grip. In contrast, the multipole ex- 
pansion technique [501 13^ is designed to model a perfect 
separation of the crack tip scale to the strip width W, i.e. 
W ^ D/vr or W ^ D/vq, respectively. In most real 
cases, crack tips are very tiny, and therefore it is theoret- 
ically desirable to describe this limit. For the phase field 
method, however, a finite strip width W is necessary, and 
a good separation of the scales therefore requires time- 
consuming large-scale calculations. We shift the system 
such that the tip remains in the horizontal center. This 
allows to study the propagation for long times until the 
crack reaches a steady state situation. Apart from this 
finite size restriction, introduced the interface width ^ as 
a numerical parameter, and the phase field method deliv- 
ers quantitative results only in the limit that all physical 
scales are much larger than this lengthscale. The latter 
has to be noticeably larger than the numerical lattice pa- 
rameter Ax, but the results show that the choice ^ = 5Aa; 
is sufficient. Therefore, to obtain quantitative agreement 
with the results from the multipole expansion method, 
we have to satisfy the hierarchy relation 

— <iy or — <W^, (76) 

vr vq 

which is numerically very hard to achieve. 
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We developed a parallel version of the phase field code 
which is running on up to 2048 processors, with sys- 
tem sizes up to 8192 x 4096 • (Ax)^ (Ax is the lattice 
unit). However, for qualitative results we typically use 
Wvr/D = 86 and D/vr^ — 1.9, where the total size of 
the system in grid points is 2048 x 800. All computations 
are performed on the supercomputer JUGENE operated 
at the Research Center Jiilich. 

The dimensionless driving force A decomposes into 
mode I and III contributions, and according to Eq. (71 1, 
it is defined for the strip geometry as 



A = A/ A/// = 



E 



51 



6hi 



2-fW \2{1 - 1^2) I 



(77) 



Here, Sj is the above mentioned fixed displacement by 
which the strip is elongated vertically, whereas Sm is a 
fixed displacement by which the strip is sheared in the z 
direction. The value A = 1 corresponds to the Griffith 
point. 



V. RESULTS 

In the following section we will give a comprehensive 
overview on the different results. As has been men- 
tioned above, we consider two different material trans- 
port mechanisms, surface diffusion (SD) and a phase 
transformation process (PT). From a theoretical point 
of view, two different physical mechanisms, i.e. viscous 
dissipation and inertial effects are important to provide 
selection mechanisms for the steady state velocity and 
the crack tip structure. These two cases can be con- 
sidered as limiting situations for slow and fast cracks, 
and here a quantitative treatment not only with phase 
field but also sharp interface methods is feasible (multi- 
pole expansion method) . The cross-over behavior, where 
both effects are relevant, is modeled using perturbation 
techniques for the multipole expansion method and fully 
dynamical phase field simulations. Furthermore, apart 
from steady state growth also branching instabilities can 
occur, which will also be discussed. Finally, we consider 
different loading modes, and we will start the discussion 
of the results for pure mode I fracture, before in the fol- 
lowing section mixed mode situations with a combination 
of mode I and mode HI loading are investigated. 

Apparently, the different physical situations and nu- 
merical methods lead to a certain complexity of the re- 
sults. A concise summary of the results is therefore given 
in Table U 



A. Opening Mode fracture 

In this section we discuss exclusively mode I fracture 
in the different variants of the model. As discussed be- 
fore, selection is of central interest for this pattern for- 
mation aspect of fracture, and two principal mechanisms 



have been introduced before, which the selection through 
viscoelastic bulk damping and the inertia limitation of 
the crack speed. In all following calculations we use 

= v = 1/3. Then, t ^ rj/E is the only remaining vis- 
cous time scale, and wc define the dimensionless viscosity 
strength x = ''-'r/'^qi where vr is the Rayleigh speed and 
vq = (L>,r-3)i/4 for the SD model, while i>o = {D/rfl'^ 
is used for the PT model. 

First, we deal with slow crack propagation with a 
steady state velocity much smaller then the Rayleigh 
speed, i.e. v <C vr. In this case dynamic effects are negli- 
gible, and the application of static elasticity is legitimate, 
X = CXI. Next, we discuss the limit of fast crack propaga- 
tion with vanishing viscous dissipation, where the steady 
state velocity v and the finite tip radius rp are selected 
by dynamic effects only, % = 0. A more general situa- 
tion, which contains both effects, will be discussed in the 
framework of a perturbation analysis using the multipole 
expansion method and fully dynamical phase field sim- 
ulations, as well as non-steady state crack growth with 
crack branching. 

The different kinetic mechanisms PT and SD lead 
to very similar results in general, apart from the fact 
that surface diffusion implies material conservation, and 
therefore the crack shapes differ (compare Figs. [2] and 
[3|. However, an important difference is that for surface 
diffusion steady state physically relevant solutions do not 
exist without viscoelastic damping. 



1. Slow cracks 

In this regime it is assumed that the sound speed is 
much larger than the crack velocity, and therefore inertial 
effects can be neglected. 

We start with reviewing the results for surface diffu- 
sion, as presented in [39]. As for all surface diffusion 
models, only multipole expansion technique results are 
available, since the modeling of surface diffusion with 
phase field methods is more cumbersome [57U59j . The 
numerical results for mode I fracture, as obtained from 
the simulations, are shown in Fig. |4] and [5] The in- 
set of Fig. [4] shows a typical steady state crack shape, 
which has the characteristic features of a finite tip radius 
and vanishing surface separation in the tail region, which 
results from the material conservation condition, as dis- 
cussed before. We point out that the crack is shown in 
the Lagrangian reference frame, i.e. without elastic dis- 
placements, as it would appear if suddenly the mode I 
loading was removed. The maximum vertical diameter 
of the crack defines here the tip scale 2h. The velocity 
plot Fig. |4] shows only finite velocities above A ss 2.6, 
and from there on a strictly monotonic increase of the 
steady state velocity as function of the driving force. We 
did not find any indications that this solution branch ter- 
minates at higher driving forces. Notice that the crack 
speed is set by the characteristic viscoelastic velocity 
vq = {DsT~^)^/'^ . Reducing the driving force coming 
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Surface Diffusion (SD) 


Phase transformation (PT) 




Viscoelastic limit: x = oo 

• Creep branch: 1 < A < 2.6 

• Regular growth for A > 2.6 

— Velocity grows monotonically with 
driving force 

- Velocity scale vo ~ (D^r"^)^/" 

- Tip scale ho ~ (Dst)^/" 

— No branching 


Viscoelastic limit: x = co 

• Creep branch: 1 < A < 2.6 

• Regular growth for A > 2.6 

— Velocity grows monotonically with 
driving force 

- Velocity scale vo ~ (D/tY^^ 

- Tip scale ho ~ (Dr)^/^ 

— No branching 


Mode I 


Inertial limit x = 

• No physically relevant solution 


Inertial limit X = 

• No self-consistently selected tip radius in 
the range 1 < A < 1.14 

• Steady state solution for A > 1.14 

— Velocity decaying function of A 

— Velocity scale vr ~ (E/p)^^^ 

— Tip scale ho ~ 

• Crack branching for A > 1.8 




Viscoelastodynamic regime < x < oo 

• Creep branch for low driving forces 

• Velocity first increases with A, then de- 
crease 

• No steady state solutions beyond a criti- 
cal driving force, afterwards branching ex- 
pected 

• Wide range of A for steady state solutions 

• Higher viscosity loads to lower crack speeds 

• Range of steady state solutions larger for 
higher viscous damping 


Viscoelastodynamic regime < x < oo 

• Creep branch for low driving forces 

• Velocity first increases with A, then de- 
crease 

• No steady state solutions beyond a criti- 
cal driving force, afterwards branching ex- 
pected 

• Small range of A for steady state solutions 

• Higher viscosity leads to lower crack speeds 

• Range of steady state solutions larger for 
higher viscous damping 


Mode III 


Viscoelastic limit x = oo 

• Strong blunting below A = 1.1 (ductile-to- 
brittle transition) 

• Steady state regime: 1.1 < A < 3.8 

— Velocity grows monotonically for 
l.K A < 3.5 

— Velocity decays for 3.5 < A < 3.8 

- Velocity scale wo ~ (-D^r"^)^/* 

- Tip scale /),„ - (D^r)^/^ 

• No steady state solutions for A > 3.8, 
where branching is expected 


ViscoelEistic limit x = oo 
• Logarithmic opening of the crack 


Mode I 
+ III 


Viscoelastic limit x = 

• Higher mode 1 contribution lead to shift of 
onset of branching towards higher A 

• Creep branch with a low velocity plateau 

• For low A faster growth for higher mode 
111 contribution may enable development of 
crack front instability 


ViscoelEistic limit x = oo 
• Logarithmic opening of the crack. 



TABLE 1: Brief summarizing comparison of the different growth modes. 
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FIG. 4: Steady state tip velocity vjvi:, for a mode I crack as a 
function of the driving force A in case of the SD model. The 
results are obtained with the multipole expansion method in 
the viscous limit. The inset shows the corresponding crack 
shape for A = 10.0. Both directions x and y are scaled with 
the half maximum height h of the crack. 
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FIG. 6: Comparison of the tip velocity v/wo as a function of 
the driving force A for mode I fracture using the PT model, 
with uo = {D/tY^'^. The solid line corresponds to the re- 
sults of the multipole expansion method with infinite viscosity 
strength. The triangles correspond to the phase field results 
with viscosity strength x = 2. The inset shows the crack 
shape for the PT model with A = 3.6 obtained with the mul- 
tipole expansion method. Both directions x and y are scaled 
with the half tail opening h of the crack. 
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FIG. 5: The crack height, as defined in the inset of Fig. |4] 
for the SD model in the viscoelastic limit. The characteristic 
lengthscale is defined as ho — (Dst)^^^. Below the value 
A = 2.6 the crack tip size becomes small, and simultaneously 
the crack velocity drops to very small values ("creep branch"); 
then almost all dissipation stems from viscous bulk damping. 



from high values, the crack velocity rapidly drops to very 
small values at A sa 2.6, and below this value the crack 
growth velocity is very close to zero (and not shown in 
the plot). Hence, there is a very rapid transition between 
different growth behaviors at this finite value of the driv- 
ing force, and we call the regime 1 < A < 2.6 the "creep 
branch" . Notice that the literal Griffith point is located 
at A = 1, but nevertheless it seems that significant crack 
growth starts only at a substantially higher driving force 
("apparent Griffith threshold"). In the creep branch al- 



most all elastic energy is dissipated by viscoelastic damp- 
ing; for a more detailed discussion of this point we refer 
to [39]. 

Fig. [5] shows the crack height h/ho with Hq = (D^t)^/'* 
as function of the driving force for surface diffusion in the 
viscoelastic limit. Again, the results are obtained by the 
multipole expansion method. The behavior is similar as 
for the crack velocity, as we also see here a monotonic in- 
crease as function of the driving force, since this increases 
the energy dissipation at the crack surfaces. When the 
driving force is reduced, the crack tip scale suddenly be- 
comes very small at A « 2.6, and below this value the 
crack becomes very sharp, h/ho <^ 1. 

Next, we discuss the results to the same regime of slow 
mode I crack growth (viscoelastic regime), but with the 
PT mechanism. Here, we performed simulations using 
both the multipole expansion approach and phase field 
modeling. Wc point out that the phase field model does 
not contain only viscoelastic damping but also inertial 
effects, i.e. the appearance of the acceleration term in 
the Newtonian equations of motion. For the phase field 
results in Fig. [H] we use x = 2, thus vr = \/2vq. This 
means that the velocity scales are not fully separated, 
and as soon as u ~ uq, the velocity has already reached a 
substantial fraction of the the sound speed, and dynami- 
cal effects start to become relevant. Since the crack speed 
is ultimately limited by the Rayleigh speed, it is therefore 
not surprising that the velocities obtained by the phase 
field method are lower than those by the multipole expan- 
sion technique, which assumes v/vji <C 1. This behavior 
is visible in Fig. [6] for driving forces A > 4. 
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Overall, the behavior for the PT model is very similar 
to the SD results, which includes in particular a mono- 
tonic increase of the steady state velocity as function of 
the driving force. Again, we did not find indications for 
crack branching at higher velocities without incrtial ef- 
fects, V ^ vji. As for the SD model, the solution branch 
with V vq terminates at A 2.6, and below we have 
again a creep branch with v <^ vq. Similarly, we also 
have here the drop of the crack tip scale h/h^ (with 
Hq = {DtY^"^ for the PT mechanism) to very small val- 
ues for A < 2.6. Here, the phase field model brings in 
another effect, which comes from the interface thickness 
^ as intrinsic numerical parameter. To obtain results that 
coincide with the multipole expansion method it is neces- 
sary to maintain the scale separation ^//i ^ 1. We have 
demonstrated for the inertial regime that it is possible to 
reach this limit, although it is numerically very demand- 
ing [38]. Here, however, we rather consider the interface 
thickness as additional "microscopic" cutoff scale, which 
prevents that the crack tip scale drops below this value; 
this effect can be seen from the steady state crack shapes, 
see Fig. [7] Therefore, for phase field modeling the creep 
branch does not exist, and consequently the crack veloc- 
ity continues to decay smoothly down to A = 1. 

For moderate driving forces around A = 4, the quali- 
tative agreement between phase field results and the ve- 
locities from the multipole expansion method is good. 

For comparison with the phase field shapes the inset of 
Fig.|6]shows a typical steady state crack contour obtained 
with the multipole expansion method in the limit of static 
viscoelasticity. Again, the shape is drawn without elastic 
displacements, which should be added to obtain the real 
shape under load. The crack tip scale is selected self- 
consistently, and the finite time cusp singularity of the 
ATG instability does not occur. Therefore, we can con- 
clude that the sole presence of viscous bulk dissipation is 
a way to cure this well-known problem |39] . Since we do 
not have a conservation law for the amount of material 
inside the crack, the pattern looks difi^erent in comparison 
to the SD crack shape. 

Finally, we remark, that within the "static" limit 
(without inertial effects) a branching instability does not 
show up for mode I fracture neither in the PT model nor 
for the SD mechanism. The phase field model of course 
contains inertial effects, and therefore for A ^ 1 one al- 
ways find branching events once v vr. In contrast, for 
pure mode HI fracture and mixed mode scenarios within 
the SD model, which will be discussed below, an insta- 
bility appears even in the static limit [39] . 



2. Inertial limit 

Here we consider situations where the crack speed be- 
comes comparable to the Rayleigh speed, while it is as- 
sumed that the viscous damping is negligible, i.e. x — 0. 
Surprisingly, for SD no physically reasonable steady state 
solutions exist, and therefore we discuss only the PT 




FIG. 7: Phase field crack shapes for x = 2 and different 
driving forces after the time t/r = 24.4; a) A = 1.25 b) 
A = 3.6 and c) A = 5.0. We set W/hg = 60.9 and ho/^ = 
2.6. The thickness of the interface corresponds the phase field 
interface width. For high driving forces, the tip radius does 
not depend on the interface thickness. Notice that for the 
lowest driving force the crack opening is not constant along 
the crack but increases towards the tail. The reason is that 
due to the elastic energy stored in the strip there is an effective 
short-range repulsion between free surfaces. 



mechanism. For that, we briefly review the results of 
our previous work 

Here, rather small scale phase field simulations '37| de- 
livered a picture, which was in conflict with very precise 
multipole expansion method results [36;, since the pre- 
dicted driving force dependence of the steady state veloc- 
ity came out with opposite slope. This discrepancy was 
resolved in [35] by performing a large series of phase field 
simulations for different system sizes W and tip scales, 
and careful extrapolation of the crack velocity to the limit 
^//i — 7- and h/W 0. The main result in this context 
was the quantitative agreement of the steady state crack 
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FIG. 8: Quantitative comparison of steady-state crack ve- 
locities obtained from the multipole expansion technique and 
the extrapolated values from phase field simulations, for the 
limiting case of x = 0. The gray line-color indicates results 
where a negative tip curvature was measured. The inset of the 
figure shows multipole expansion crack shapes of the stable 
(black curve A — 1.3) and the unstable solution (gray curve 
A — 2.3). Both directions x and y are scaled with the half tail 
opening h of the crack. Below the point Ac ~ 1.14, indicated 
by the dotted line, we show the velocity of the dissipation-free 
solution, where the tip radius ro is selected by a microscopic 
length scale. 



velocities v/vji obtained from both numerical methods, 
as shown in Fig. [8| The small deviation for A = 1.8 
is due to the fact that this value is already close to the 
threshold of the branching instability, which cannot be 
captured by the multipole expansion method. With this 
costly quantitative comparison, we found in particular 
evidence for the remarkable prediction that the steady 
state velocity decays weakly with increasing driving force 
|36j . Nevertheless, the product vh/D, which controls the 
dissipation, is still growing monotonically. This counter- 
intuitive outcome means that within the dynamic limit 
ix — 0) of the model the dissipation is mainly increased 
due to tip blunting instead of a rise of the crack speed. 
Tip blunting then always leads to a tip branching in- 
stability for higher driving forces, due to a secondary 



ATG-instability as mentioned in section III In the multi- 



pole expansion method, which captures only steady state 
solutions, this transition towards unstable crack growth 
is reflected by a change of sign of the tip curvature at 
A « 1.8, which is in agreement with the critical driving 
force for the branching instability in phase field model- 
ing. In Fig. |8]we indicate this change by a change from 
the black to the gray line-color, and in the inset we show 
two corresponding crack shapes. For further details we 
refer to |36j . 

Although the model provides a selection of the crack 
tip scale and velocity in the limit of vanishing dissipa- 
tion, it suffers from the fact that for the small range of 



driving forces near the Griffith point the velocity of the 
crack is finite while the size of the crack tip approaches 
zero. More precisely, the velocity branch in Fig. |8] termi- 
nates at A — 1.14 with a finite velocity, and below this 
value no steady state solutions are found. Thus, the so- 
lution branch does not naturally connect to the Griffith 
point A = 1. Here, another tip scale mechanism would 
be necessary, in order to restore selection. In the phase 
field method (without performing the extrapolation to 
the sharp interface limit £^/h — 0) the interface thick- 
ness serves here as numerical tip scale selection mecha- 
nism |37| . In the same way, setting the minimal allowed 
opening by hand, also the situation for the multipole ex- 
pansion method can be improved [36] . However, since our 
intention was to formulate a continuum model of fracture 
which is independent of any microscopic length scales, for 
the present purpose, such introduction of a finite cutoff 
length scale is unsatisfactory. Then, on the other hand, 
the sudden velocity jump and the subsequent decay with 
increasing driving force are unavoidable outcomes and 
seem to disagree with intuitive expectations. This, to- 
gether with the fact, that in the inertial limit no steady 
state solutions exist for SD, has been a major motiva- 
tion for considering viscoelastic bulk dissipation as an 
alternative selection mechanism. We point out that the 
selection of tip velocity and radius via viscous bulk dis- 
sipation as discussed in the previous subsection neither 
suffers from the problem of finite velocities slightly above 
the Griffith point nor requires the introduction of an ad- 
ditional microscopic cutoff length scale. On the other 
hand, for mode I the tip branching instability does not 
occur without inertial effects. Therefore, to describe the 
full picture of crack propagation under mode I loading, 
it is highly desirable to account for both viscous dissipa- 
tion as well as dynamic effects, as will be discussed in the 
next subsection. 



3. Viscoelastodynamic regime 

The rigorous treatment of the regime where both in- 
ertial and viscoelastic effects are relevant has been per- 
formed with the phase field model, and additionally an 
approximative description with the multipole expansion 
method is possible. The model, which has been intro- 
duced in section |IV A3[ has the advantage that it is 
exact without viscous damping in the inertial regime, 
X = 0, and it is a rigorous perturbation expansion for 
slow cracks, w <C w^, using the viscosity vr/h as small 
expansion parameter, where h is the crack tip scale, hence 
operating in the regime of large values of x — > oo. There- 
fore, this model allows to gain qualitative insights into 
the full problem of dynamic mode I crack propagation 
including viscous bulk friction, both for the PT mecha- 
nism and SD. 

We start the discussion of the results again for the SD 
model, and the results are all obtained by the multipole 
expansion method. As mentioned above, no physically 
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FIG. 9: Viscoelastodynamic model results for the steady state 
velocity v/vr as a function of the driving force A in case of 
mode I crack propagation using the SD model. The solid line 
corresponds to x = 1-0, the dashed one to x = 0.5. 



reasonable steady state solutions exist in the purely iner- 
tial limit, x = 0; we will return to this point below. 
However, with the inclusion of viscous effects within the 
present model, steady state solutions exist for finite val- 
ues of Xj as shown in Fig. [9j In this plot the velocity is 
shown on the scale of the Rayleigh speed vr^ which equals 
the viscous scale uq for x = 1- As in the purely viscous 
limit, fast growth does not start at the literal Griffith 
point A = 1, but a higher value, which is located at 
around A w 2.1 for both shown paremeters x = 0.5 and 
X ~ 1 (below this value we have again the creep branch 
with V <^ vq). From this point on the velocity first in- 
creases with increasing driving force until it reaches a 
maximum and then it again decreases until the termina- 
tion in a bifurcation. We expect crack branching beyond 
this point, since steady state solutions do not exist in this 
regime. 

In agreement with the intuitive expectation that the 
cracks should slow down by viscous damping, the velocity 
decreases with the increase of the viscosity strength x = 
vj^/^/DgT^-^. At the same time, the bifurcation point, 
beyond which no steady state solutions exist, is shifted 
to higher driving forces. This fact indicates that viscous 
damping suppresses crack branching, in agreement with 
the earlier observation that it does not occur in the purely 
viscous limit. For very strong viscous damping, x ^ oo, 
the velocity becomes small and is then set by the viscous 
speed scale vq, as depicted in Fig.|4j where it is a purely 
monotonically growing function of the driving force. 

On the other hand, the results show that upon reduc- 
tion of the viscous strength, i.e. for smaller values of x 
the velocity increases. Finally, the curves first touch and 
then terminate at the Rayleigh velocity, which is the up- 
per theoretical limit for mode I fracture. Then, for the 
inertia limit the curves would start with the decaying 
part of the curve at a finite value of A > 1 but with 
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FIG. 10: Viscoelastodynamic model results for the steady 
state tip velocity v/va as a, function of the driving force A in 
case of mode I crack propagation using the PT model. The 
solid black line corresponds to x = 1-0, the dashed to x = 0.5. 
The gray colored solid line shows the steady state velocities 
in the inertial limit, when x = 0.0 (see Fig. [8|. 



V — Wfl, and we do not consider this as a physically plau- 
sible solution. However, in the framework of the present 
model it becomes thereby understandable why no "rea- 
sonable" solutions exist in the elastodynamic limit. 

Next, we inspect the behavior of the model for PT 
dynamics, as obtained from the multipole expansion 
method. The results for the same two different values 
of X are shown in Fig. 



10 and exhibit a qualitatively 



similar behavior as for SD. Also here, the growth starts 
from an "apparent" Griffith point, which coincides with 
the value for SD, and below we find "creep solutions" 
with very low velocities, which are not shown in the plot. 
The reason for the agreement of these apparent Griffith 
points is that for v <^ vq the chemical potential is ba- 
sically zero along the crack, both for SD and PT; also, 
below this value of the driving force the behavior is dic- 
tated by bulk dissipation, and therefore the transport 
mechanism on the crack surfaces plays only a minor role. 

From the apparent Griffith point on the velocities in- 
crease monotonically quite rapidly up to a maximal value. 
Then the velocity maximum is followed by a small range 
of driving forces, where the velocity decreases with in- 
creasing driving force. Crack branching is expected be- 
yond the bifurcation point. With increasing viscosity 
strength x the driving force of maximal velocity as well 
as the point where the curvature turns negative are both 
shifted to higher driving forces. This supports the con- 
clusion that dynamic effects favor of the occurrence the 
tip splitting instability, or - vice versa - the presence of 
viscous bulk friction helps to stabilize the crack against 
the tip splitting instability, which is also qualitatively 
supported by fully dynamic phase field simulations. 

For X — oo the results come closer to the previous 
curve for the viscoelastic limit of the PT model, as shown 
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FIG. 11: Irregular tip splitting scenario for a viscosity 
strength of x = 2 and A = 10.0. We set Wvr/D = 170 
and D/vr^ — 1.9, and the size of the system is 1600 x 4096 
grid points. The time t is given in units D/v\. The thickness 
of the interface corresponds the phase field interface width. 
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FIG. 12: Steady state propagation velocity as function of the 
driving force for pure mode III and a mixture with Aj/j /A = 
0.85 are displayed, for the SD model. The gray line belongs to 
steady state solutions with negative tip curvature. The inset 
of the figure shows the steady state crack shapes of the stable 
and the unstable solution in the case of mixed mode loading 
with A//j/A = 0.85 and a total driving force of A = 3.6. 
Both directions x and y are scaled with the half maximum 
height h of the crack. 



B. Mixed Mode fracture 



in Fig. [6] where the velocity is then on the scale Wq. In 
contrast, with a decrease of the viscosity strength the 
point of maximal velocity is shifted more and more to 
the left until, finally, in the case of vanishing viscosity, 
X = 0, only the decaying part remains (see also Fig. [8|. 

The PT and SD model have in common that the crack 
velocity decreases with the increase of viscosity strength 
X, while the turning point is shifted to higher driving 
forces. A difference is, that the point of maximal velocity 
and moreover the tip curvature turning point appear at 
much higher driving forces for the SD than for the PT 
model, which especially also enlarges the regime, where 
the velocity decreases with increasing driving force. 

The phase field method allows to study crack growth 
also in regimes, where no steady state solutions exist. 
However, a quantitative determination of the onset of the 
branching instability is computationally very expensive. 
The onset of the irregular tip splitting behavior depends, 
in particular, sensitively on the system size, because in 
relatively small systems the branches of the crack can- 
not separate since they are repelled by the boundaries. 
Therefore, the steady state growth is always stabilized by 
finite size effects. On the other hand, initial conditions 
can trigger an instability, and then a long transient is 
required to get back to steady state solutions. However, 
as shown in Fig. |1H even for a relatively high viscosity 
strength x = 2 we find the irregular tip splitting behavior 
for A = 10.0. 



Here we discuss predictions of the model beyond a pure 
mode I loading. It turns out that the results change 
even qualitatively for mode HI. For a broader picture, 
we study also situations with mixed mode loading (mode 
I + mode HI), but still assume that the crack shape is 
translational invariant in the direction of the crack front 
line. We therefore suppress effects like the development 
of a helical instability as studied in [6T. We limit our 
investigations here to viscous effects (i.e. v vr), and 
only for the phase field model we also take into account 
inertial effects. As before, we make the simplifying as- 
sumption = C = 1/3, and therefore viscosity introduces 
only one additional time scale, r = t]/E. 

First, we review our findings concerning the crack be- 
havior in the mixed mode scenario for the case of the 
SD model [33], which were obtained by the multipole ex- 
pansion method. As shown in Fig. [T2j the crack speed 
increases with the driving force for pure mode III, until 
it reaches a maximum at A « 3.5, then it decreases, and 
obviously steady state solutions do not exist beyond the 
point A « 3.8, where the stable branch merges with an- 
other (unstable) solution. Beyond the bifurcation point 
A w 3.8 we expect crack branching, in analogy to our 
findings for fast dynamic mode I fracture, as discussed 
in the previous section. It is quite remarkable, that the 
presence of mode HI loading contribution leads to the oc- 
currence of the tip branching instability even in the case 
of static elasticity. 



In Fig. 13 we also show the maximum height of the 



21 



1.5 



0.5 




pure mode III 
85% mode III 



1 2 3 4 5 6 

A 

FIG. 13: Half crack height as function of the driving force for 
pure mode III and a mixture with Aj_r//A — 0.85, for the 
SD model. The lengthscale used here is ho = (Dr)^''*. The 
gray line belongs to steady state solutions with negative tip 
curvature. 



crack as function of the driving force for different load- 
ings. At A « 1.1 the size of the mode III steady state 
crack diverges and u — > 0. The viscous dissipation be- 
comes negligible here, but the surface dissipation remains 
finite. This point can be interpreted as the point of 
ductile-to-brittle transition: Below it the size grovi^s in- 
definitely in time and the crack slows down, while above 
this point steady state solutions with a finite tip scale 
exist. 

Starting from a pure mode III crack, we can now in- 
clude additional mode I loadings. Fig. 12 shows that this 



shifts the bifurcation point towards higher values and 
therefore extends the range of steady state solutions to- 
wards higher driving forces. From this we can conclude 
that mode III contributions favor the appearance of the 
tip splitting instability. In contrast, the preceding re- 
sults suggets that inertial effects should push the onset 
of branching back towards lower driving forces. 

It is important to note that mode I and mode III have 
a different behavior, which is due to the behavior of the 
stresses on the crack surfaces. We focus here on the 
elastic fields far away from the tip, and in this region 
the behavior is purely elastic. Without inertial effects, 
normal and shear stresses vanish on the crack surfaces, 
and therefore it is the tangential stress component which 
determines the elastic contribution to the chemical po- 
tential. From the singular contributions to the stress 
[H [TS] one obtains on the crack lips elastic chemical po- 
tential contributions iiei{x — ^ — oo) ~ l/x for mode III 
and Hei ~ l/x"^ for mode I. This weaker decay of the sin- 
gular fields for mode III also influences the crack shapes. 

Let us look at the asymptotic shape y{x) of a crack in 
the SD model (in the tail region x — >■ — oo), and focus on 
polynomial terms. As discussed in |35j exponential terms 
are very important for the selection of the crack velocity 
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FIG. 14: Qualitative comparison of crack tip velocities v/vq 
as a function of driving force A for 50% mode III loading 
in the case of the PT model. The solid line corresponds to 
the results of the multipole expansion method in the viscous 
limit. The symbols correspond to the phase field results with 
a viscosity strength x = 2. For the phase field, we used 
Wvr/D — 86 and D /vr,£^ — 1.9; the size of the system in grid 
points is 2048 x 800. The inset shows the multipole expansion 
method steady state crack shape for a total driving force of 
A = 12.0. Both directions x and y are scaled with the half 
tail-opening h of the crack. 



and tip scale, which is related to the supression of growing 
exponentials in the tail region. Remaining exponentially 
decaying terms are small in comparison to the power law 
terms in the asymptotic regime, and therefore we sup- 
press them here. By the steady state equations of mo- 
tion Eqs. Q and ^ we obtain vy{x) ~ Dgy'" — const /x"^ 
for mode III. Therefore, we obtain a scaling behavior for 



surface diffusion as y{x) 



Correspondingly, in the 



case of pure mode I loading for SD the shape function 
decays hke y{x) ^ x^^, which is substantially faster. 

In the case of the PT model this effect is more pro- 
nounced. Since the amount of "material" inside the 
crack is not conserved, it is reasonable to look for shape 
functions y(x) = h + Sy{x), again with purely polyno- 
mial functions Sy{x). Neglecting the second derivative 
dy"{x) from the curvature contribution to the chemical 
potential, we obtain from Eqs. (10 1 and (11) in the co- 



moving frame of reference the following ordinary differen- 
tial equation in the asymptotic regime: —v6y'{x) ~ 1/x. 
Hence, in the case of a finite mode III contribution the 
shape function does not even decay but instead weakly 
grows like y{x) ~ ln(x). This slow opening of the crack 
becomes negligible for higher crack speeds. This weak 
logarithmic growth of the asymptotic crack shape is also 
confirmed by the multipole expansion method simula- 



tions as shown in the inset of Fig. 14 in the case of 50% 
mode III contribution. 

For the phase field simulations, which are also shown 
in Fig. 14 for A// A/// = 1 and x = 2, we find a re- 
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markably good agreement with the sharp interface re- 
sults over the whole range of driving forces. With the 
phase field method, we observe two different kinds of 
growth: Slightly above the Griffith point up to a driv- 
ing force of about A « 1.4, we obtain solutions with al- 
most zero velocity and an asymptotically growing crack 
opening similar to what the shape from the multipole 
expansion method shows. Then above this point the so- 
lutions seem to regularize, and this weak growth of the 
shape function is no longer observable using the phase 
field method, probably due to the higher growth veloc- 
ity. For pure mode III loading this transition point is 
shifted to an even higher driving force of about A « 2.0, 
as we analyzed by means of phase field simulations. 

We point out that these shape peculiarities can also 
be interpreted from the more general argument that no 
stationary shapes exist in mode III if only elastic effects 
are taken into account. For this loading mode the only 
two nonvanishing stress components a^z and can be 
expressed as real and imaginary part of an analytic func- 
tion. An equilibrium solution would require a vanishing 
chemical potential along the entire crack shape, which in 
turn demands that the aforementioned analytic function 
is zero there, since the elastic part is quadratic in both 
stress components. If an analytical function is zero along 
a line segment, it must vanish everywhere; this, however, 
is not compatible with a nontrivial remote stress field. 



VI. SUMMARY AND CONCLUSION 

We have presented a continuum description of fracture 
in the spirit of elastically driven interfacial pattern forma- 
tion processes. This description leads to moving bound- 
ary problems, where not only the propagation velocity 
but also the entire shape, and especially the tip radius, 
have to be selected self-consistently. In particular, we 
have discussed two different mechanisms of crack propa- 
gation: In the first case the crack is considered to advance 
by material diffusion along the crack surface. Secondly, 
we interpret fracture as a first order phase transformation 
process of the solid material to a "dense gas" phase. 

Scaling arguments were given that - to cure the fi- 
nite time cusp singularity of the ATG instability - one 
necessarily needs an independent selection of the tip ra- 
dius and the steady state velocity, which is not possi- 
ble if solely static linear elasticity is taken into account. 
Therefore, apart from capillarity and linear elasticity, ad- 
ditional physical effects are required for the determina- 
tion of additional length and time scales. Since we focus 
on gaining fundamental insights into the phenomenon of 
fracture, we have concentrated here only on well estab- 
lished theoretical concepts for dynamic elasticity and vis- 
cous dissipation. 

The arising moving boundary problems have been 
solved by two complementary methods. First, we 
have developed an efhcient steady state sharp interface 
method based on the expansion of the mechanical state 



in eigenfunctions of a straight mathematical cut. Second, 
a fully dynamic phase field description of crack propaga- 
tion by first order phase transformation processes has 
been developed. 

With these numerical tools at hand, we have obtained 
profound insights into the model behavior of our con- 
tinuum description of fracture. In particular, we have 
extensively discussed mode I fracture, where the coupled 
influence of dynamic elasticity and viscous dissipation 
leads to a model behavior which reproduces three im- 
portant generic features of fracture: (i) The saturation 
of steady state velocities appreciably below the Rayleigh 
speed, (ii) parameter regimes, where the steady state ve- 
locity increases with increasing driving force, and (iii) 
a crack branching instability for high applied tensions. 
Apart from this also mixtures of mode I and mode III 
loadings have been discussed, and we have found in par- 
ticular a different behavior of the crack shapes, as well 
as a change of the branching behavior. The main results 
are summarized in Table U 

The future challenge is to combine the ideas, ap- 
proaches and results from this work to "conventional" 
fracture models, to address the question of the relevance 
of the different physical mechanisms. From a thermody- 
namical point of view, there is a driving force for elasti- 
cally induced phase transformations, which leads to the 
fracture models presented here. However, it has to be 
addressed to which extend, and under which environ- 
mental conditions, these processes can be understood as 
dominant mechanisms for crack propagation, i.e. for in- 
stance the diffusion along crack surfaces can be an effi- 
cient mechanism in comparison to the pure bond break- 
ing. We note, that due to the small tip scales the mate- 
rial transport is necessary only on very short distances, 
and therefore a mechanism like surface diffusion, which 
is usually assumed to be slow, can still lead to fast crack 
propagation. Here it should be pointed out that on such 
small scales a pure continuum description may not be 
quantitatively accurate, but can still capture the essen- 
tial physical mechanisms. Furthermore, recent experi- 
mental investigations of fracture in brittle gels possibly 
reveal macroscopic scales |62| . 

In general, the question concerning energy barriers 
should play a central role and should shed light on the 
relevance of the different mechanisms. We expect that 
material transport should become relevant at elevated 
temperatures. In the conventional picture for brittle ma- 
terials a few bonds per atom have to be broken to advance 
the crack by one lattice unit, and this event takes place 
very localized at the (sharp) tip. The energetic cost for 
such an "event" is on the order of eV/atom. In contrast, 
for the "material transport picture" the overall energetic 
expense is the same (since the same amount of new inter- 
faces is created) , but a change of the bonding situation 
is required for several atoms. However, since the diffus- 
ing atoms do not have to be completely detached from 
the surfaces, energetically efficient low-barrier paths may 
exist for the motion to the next lattice site, and there- 
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fore the effective diffusion constant can be relatively high. 
Furthermore, surface reactions, as a recently predicted 
amorphization of diamond, can lead to a bond weaken- 
ing and material softening of even very brittle materials 
within short times [631 and could facilitate even higher 
transport rates. 

We notice that in some cases nonlinear elastic correc- 
tions may play an important role |481I49| and even lead to 
a high-speed oscillatory instability [SU] ■ For more ductile 
materials, plastic processes due to dislocation emission 
are important, and they have not yet been taken into 
account. We expect bulk dissipation through plasticity 
to play a similar role as viscous damping, as has been 
demonstrated above. Apart from that, these theories in- 
troduce the concept of a yield stress ay , which is a natural 
cutoff for the stress singularity. Therefore, from point of 
view of crack tip selection, one can expect the radius tq 
to be determined by this cutoff, i.e. tq ~ K^/ay. Since 
this leads to strong blunting, surface diffusion is probably 
not efficient for high driving forces, and in fact the con- 



tribution of surface diffusion to the propagation velocity 
would be a decaying function of the driving force, hence 
only for small A it may compete with a bond-breaking 
mechanism together with plastic flow. Another aspect is 
that plastic effects lead to large deformations, and there- 
fore from a technical point of view it would then be desir- 
able to describe the material transport processes in the 
deformed system, which suggests the use of a Eulerian 
rather than a Lagrangian description. 
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